/*
 * test.cpp
 *
 *  Created on: Mar 16, 2012
 *      Author: pickrell
 */

#include "ManyTrees.h"
using namespace std;

int main(){

	//ManyTrees t("sim_nomig_test2.trees.gz");

    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_ranlxs2;
    r = gsl_rng_alloc(T);
    int seed = (int) time(0);
    //seed = 103484;
    gsl_rng_set(r, seed);
    //Alignment t("test_3haps.gz");
	//Alignment t("/Users/pickrell/projects/admix_sims/real_data_1kg/vcf_all_20101123/chr9_95092568_95161740.tree_in.gz");
    //t.compress();
	//t.print("testaln.gz");

	PhyloPop_params p;
	//Alignment a("test_4haps.gz");
	Alignment a("treein40.gz");

 	TreeState ts(&a, &p);
 	double lk = ts.tree->set_random_UPGMA(&a, 20000, 20000.0, 1e-8, r);
 	cout << ts.tree->get_newick_format() << "\n";
 	cout << "lk "<< lk <<"\n";
 	ts.tree->print("test40");
 	//ts.set_to_newick("chr9_95092568_95161740.treeout.gz");
 	//ts.init_tip_labels(0.0001, r);
 	//t.tree->print("testout");
 	//ManyTrees t("popsplit_in.tmp", &p);
 	//ManyTrees t("testin4", &p);
	//t.set_bpoints(50);
	//t.bpoints[38] = 0.18;
 	cout << ts.tree->get_loggamma(0.0354042, 3.33333,  0.010101) << "\n";//2.99477
 	cout << ts.tree->get_loggamma(1.83115e-317 ,0.00301114 ,0.25) << "\n"; //inf?
 	cout << ts.tree->get_loggamma(0, 0.00223964 ,0.333333) << "\n"; //-inf?
 	/*
 	t.set_bpoints(2);
	t.bpoints[0] = 0.18;
	for (int i = 0; i < t.bpoints.size(); i++){
		cout << i << " "<< t.bpoints[i] << "\n";
	}
 	t.init_popsizes();
 	t.optim_all_popsizes();

 	vector<double> tmp;
 	for (int i = 0; i < t.popsizes.size(); i++){
 		tmp.push_back(t.popsizes[i]);
 	}
 	t.print_state("testout");
 	t.init_structure(0, r);
 	t.set_known(51);
 	for (int i = 0; i < t.popsizes1.size(); i++){
 		t.popsizes1[i] = 0.5;
 		t.popsizes2[i] = 2.5;
 	}
	for (int i = 0; i < t.popsizes3.size(); i++){
 		t.popsizes3[i] = 0.5;
 	}
 	t.print_all_trees("test");
 	//t.optim_all_struct_popsizes();
 	cout << "llk1 "<< t.d[0].structured_llk_no3(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2) <<"\n";

 	vector<Graph::vertex_descriptor> pos = t.d.at(0).tree->get_postorder_traversal(t.d.at(0).tree->root, 100);
 	t.d.at(0).move_label_subset(pos, 1);

 	for (int i = 0; i < t.popsizes1.size(); i++){
 		t.popsizes1[i] = tmp[i];
 		t.popsizes2[i] = tmp[i];
 	}
 	t.p1 = 0.2;
 	t.print_all_trees("test1");
 	cout << "llk2 "<< t.d[0].structured_llk_no3(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2) <<"\n";
*/









	//cout << t.llk_structured_no3() << "\n";

 	//t.optimize_all_tips();
 	//cout << t.llk_structured_no3() << "\n";

 	//t.optim_all_struct();
 	//t.optim_all_struct();
 	//t.print_struct_state("testout.struct");
 	//t.print_all_trees("test1");

 	//TRUTH
/*
 	for (int i = 0; i < t.popsizes1.size(); i++){
 		t.popsizes1[i] = 0.25;
 		t.popsizes2[i] = 0.25;
 	}
 	for (int i = 0; i < t.popsizes3.size(); i++){
 		t.popsizes3[i] = 0.5;
 	}
 	t.popsizes2[t.popsizes2.size()-1] = 2.5;
 	t.popsizes2[t.popsizes2.size()-2] = 2.5;
*/
 	//t.optim_all_struct();
 	//t.d[2].move_label_subset(inorder, 2);
 	//cout << "llk2 "<< t.d[2].structured_llk_no3(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2) <<"\n";

	//;t.d[2].tree->print("before");
 	//t.optimize_all_tips();
 	//t.optimize_all_tips();

	//pair<Graph::vertex_descriptor, Graph::vertex_descriptor> ch = t.d[2].tree->get_child_nodes(t.d[2].tree->root);
	//Graph::vertex_descriptor c1 = ch.first;
	//if (t.d[2].tree->g[c1].pop != 3) c1 = ch.second;
	//ch = t.d[2].tree->get_child_nodes(c1);
	//c1 = ch.first;
	//if (t.d[2].tree->g[c1].pop != 3) c1 = ch.second;
	//t.d[2].optimize_single_tip(c1, t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2);
	//cout << "llk2 "<< t.d[2].structured_llk_no3(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2) <<"\n";
	//t.d[2].tree->print("after2");

	//for (int i = 0; i < t.popsizes1.size(); i++){
 	//	t.popsizes1[i] = 0.5;
 	//	t.popsizes2[i] = 0.5;
 	//}
 	//for (int i = 0; i < t.popsizes3.size(); i++){
 	//	t.popsizes3[i] = 0.5;
 	//}
 	//t.popsizes2[t.popsizes2.size()-1] = 5;
 	//t.popsizes2[t.popsizes2.size()-2] = 5;
 	//t.p1 = 0.5;
 	//t.optimize_all_tips();
 	//t.optimize_all_tips();
 	//t.d[2].tree->print("after3");
 	//cout << "llk3 "<< t.d[2].structured_llk_no3(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2) <<"\n";
 	//t.d[2].iterate_optimize_all_tips(t.p1, &t.bpoints1, &t.popsizes1, &t.popsizes2);


 	//t.d[0].init_tip_labels(0.25, r);
 	//t.d[0].tree->print("before");
 	//vector<double> bp1;
 	//bp1.push_back(0.05); bp1.push_back(0.15); bp1.push_back(0.25);
 	//vector<double> p1;
 	//vector<double> p2;
 	//p1.push_back(0.5); p1.push_back(0.5); p1.push_back(0.5);
 	//p2.push_back(0.5); p2.push_back(0.5); p2.push_back(0.5);


 	//t.p1 = 0.5;
 	//cout << t.d[0].structured_llk_no3(t.p1, &bp1, &p1, &p2) << "\n";
 	//t.d[0].n1 = 3;
 	//t.d[0].n2 = 20;
 	//cout << t.d[0].structured_llk_binom(0.2) <<"\n";
 	//t.set_bpoints(50);
 	//t.init_popsizes();
 	//t.optim_all_popsizes();
 	//t.init_structure(48, r);
 	//t.d[0].print_coals12();

 	//t.d[0].optimize_all_tips(t.p1, &t.bpoints1,  &t.popsizes1, &t.popsizes2);
 	//t.optimize_all_tips();
 	//t.print_all_trees("test");
 	//t.d.at(0).tree->print("tree1");
	//t.d.at(1).tree->print("tree2");
 	//cout << t.llk_structured() << "\n"; cout.flush();
 	//t.iterate_optim_all_struct();
 	//cout << t.llk_structured() << "\n"; cout.flush();
 	//t.print_struct_state("teststruct");
 	//t.optim_all_struct();
 	//cout << t.llk_structured() << "\n"; cout.flush();
 	//cout << "here1\n"; cout.flush();
 	//t.d.at(0).tree->print("before");
 	//t.d.at(0).optimize_all_tips(&t.bpoints1, &t.bpoints2, &t.popsizes1, &t.popsizes2);
 	//t.d.at(1).optimize_all_tips(&t.bpoints1, &t.bpoints2, &t.popsizes1, &t.popsizes2);
 	//cout << "here\n"; cout.flush();


	//t.d.at(0).iterate_optimize_all_tips(&t.bpoints1, &t.bpoints2, &t.popsizes1, &t.popsizes2);
 	//cout << "here\n"; cout.flush();
	//cout << t.llk_structured() << "\n";
	//t.d.at(0).tree->print("after2");
    //ts.add_pop(); exit(0);
    //ts.tree->print("before");
    //ts.many_local_hillclimb();
    //ts.tree->print("after");
    //vector<bool> t2 = t.get_hap("h2");
    //vector<string> allhapnames = t.list_haps();
	//vector<string> startpops;
	//startpops.push_back(allhapnames[0]); startpops.push_back(allhapnames[1]); startpops.push_back(allhapnames[2]);
	//vector<vector<bool> > starthaps;
	//starthaps.push_back(t.get_hap(startpops[0])); starthaps.push_back(t.get_hap(startpops[1])); starthaps.push_back(t.get_hap(startpops[2]));
	//PopTree tt(startpops, starthaps, t.len, &p);

	//vector<Graph::vertex_descriptor> postorder = ts.tree->get_postorder_traversal(3);

	//for (vector<Graph::vertex_descriptor>::iterator it = postorder.begin(); it != postorder.end(); it++) cout << tt.g[*it].index << " ";
	//cout << "\n";
	//tt.compute_tree_llk();
	//ts.tree->print("before");
 	//vector<bool> newhap = t.get_hap("hap31");
 	//Graph::vertex_descriptor vadd = ts.tree->add_tip(postorder[1], newhap, "hap31");
	//ts.add_pop();
	//ts.tree->print("before");
	//postorder = ts.tree->get_postorder_traversal(4);
	//ts.tree->local_rearrange(postorder[4], 2);
	//ts.tree->print("after");
	//ts.iterate_hillclimb();
	//ts.tree->compute_tree_llk();
 	//cout << ts.tree->current_llk << "\n";


	//
	//
	//ts.tree->print("after");
	//ts.tree->compute_tree_llk(vadd);
	///cout << ts.tree->current_llk << "\n";
	//Graph::vertex_descriptor v = postorder[3];
	//tt.g[v].height = 0.45;
	//tt.set_ls(v);
	//tt.set_ls(postorder[3]);
	//tt.compute_tree_llk();
	//cout << tt.current_llk <<"\n";
	//for (int i = 0 ; i < 100; i++){
	//	double t = (double)i/100.0;
	//	double t2 = 0.02*t;
	//	tt.g[v].height = t2;
	//	tt.compute_tree_llk();
	//	cout << i << " "<< t2 << " "<<  tt.current_llk << "\n";
	//}
	//tt.optimize_height(postorder[postorder.size()-1]);
	//tt.optimize_height_gs(postorder[3]);
	//tt.optimize_height_gs(postorder[postorder.size()-1]);
	//tt.optimize_height(postorder[3]);
	//vector<Graph::vertex_descriptor> tv;
	//tv.push_back(postorder[3]);
	//tv.push_back(postorder[postorder.size()-1]);
	//tt.optimize_heights(tv);
	//tt.print("before");
	//tt.compute_tree_llk();
	//cout << tt.current_llk << "\n";
	//cout << tt.ntip << " ntip\n";
	//tt.add_tip(postorder[3], t.get_hap("h4"), "h4");
	//cout << tt.ntip << " ntip\n";
	//tt.local_rearrange(postorder[3], 2);
	//tt.print("before");
	//tt.move_root(4);
	//tt.optimize_all_heights();
	//cout << tt.current_llk << "\n";
	//tt.print("after");
	//tt.compute_tree_llk();
	//cout << tt.current_llk << "\n";
	//tt.compute_tree_llk(postorder[3]);
	//cout << tt.current_llk << "\n";

	return 0;
}

